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A general lattice Boltzmann method for simulation of fluids with tailored transport coefficients 
is presented. It is based on the recently introduced quasi-equilibrium kinetic models, and a general 
lattice Boltzmann implementation is developed. Lattice Boltzmann models for isothermal binary 
mixtures with a given Schmidt number, and for a weakly compressible flow with a given Prandtl 
number are derived and validated. 
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I. INTRODUCTION 

The lattice Boltzmann (LB) method is a powerful approach to hydrodynamics, with applications ranging from 
large Reynolds number flows to flows at a micron scale, porous media and multi-phase flows [l|. However, in spite of 
its rapid development over a decade, the method is still lacking a systematic and flexible construction of numerical 
schemes for situations beyond a simple fluid, such as mixture models with a given Schmidt number, or thermal models 
with a given Prandtl number. 

In this paper, we introduce a general method of constructing the lattice Boltzmann models. Our approach is 
based on recently introduced quasi-equilibrium (QE) kinetic models @, The structure of the paper is as follows: 
First, for the sake of completeness, we remind the construction of continuous time-space QE models. After that, we 
derive a lattice Boltzmann discretization scheme of the QE models. The resulting quasi-equilibrium lattice Boltzmann 
method (QELBM) is illustrated with two examples of particular interest. The first example is the thermal model 
with a prescribed Prandtl numbers. The second example is the isothermal binary mixture model with a prescribed 
Schmidt number. Basic steps of construction of these models are presented and numerical validation is provided. 

II. QUASI-EQUILIBRIUM KINETIC MODELS 
A. General 

We address the class of lattice Boltzmann models equipped with the Boltzmann entropy function of the form, 

where f% > are populations of the discrete velocities Vi, i = 1, . . . , n, and Wi > are correspondingweights. A wide 
class of relevant entropy functions ([I]) pertinent to simulation of hydrodynamics was described in H, EE 0,1- In 
the examples below we shall use the 7J-functions for the isothermal models Q and the recently introduced weakly 
compressible thermal models @ in two dimensions, for which the set of nine velocities Vi and corresponding weights 
Wi are 

v x = {0,1,0,-1,0,1,-1,-1,1} 

v y = {0,0,1,0,-1,1,1,-1,-1} (2) 
W = (1/36) {16, 4, 4, 4, 4, 1,1, 1,1}. 

We begin with a generic construction of the quasi-equilibrium kinetic models 2, 3] specified for the discrete velocity 
case. Let M = {Mi , . . . , Mk M } be a set of locally conserved fields, and N = { Ni , . . . , Nk N } be a set of quasi-conserved 
slow fields. Functions M m and Ni are assumed linear functions (moments) of the populations, M m = J^. M m ifi and 
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Ni = J2i Nufi- The choice of M and N depends on the particular problem but we always assume that the density is 
included into the list M, p = Mi = Yl7=i The equilibrium population vector f cq (M) is defined as the minimizcr 
of the i?-function ([lj under fixed M. The quasi-equilibrium f*(M,N) is defined as the minimizer of H under fixed 
M and N. By construction, functions /* and / cq satisfy consistency relations 

M(f*(M,N)) = M, 

N(f*(M,N)) = N, 

M(f cq (M)) = M. (3) 

The quasi-equilibrium kinetic model reads 



d f M + v ia d a fi = - [.n f*(M(f),N(f))} - [f*(M(f),N(f)) fP(M(f))} . (4) 
n t 2 

Denoting the right hand side (the collision integral) as Qi, it is easy to see that consistency condition implies local 
conservation laws, M (Q) — 0. Note that the first part of the collision integral which describes relaxation to the quasi- 
equilibrium, also conserves the iV-fields, N(f — /*) = 0. Moreover, it is straightforward to prove the i?-theorem. For 
that, it suffices to rewrite 

Qi - ~ [fi f?] - ( [Si - m ■ (5) 
The entropy production a = 'dfi)Qi becomes 

'~4t-(£KM^)g ta (& 



Ui m, (6) 



and is non-positive semi-definite provided relaxation times satisfy the hierarchy, 

n < t 2 . (7) 

Thus, in the QE models, relaxation to the equilibrium is split in two steps. In the first step, the distribution functions 
relaxes to the quasi-equilibrium with the faster relaxation time t\. In the second step, the quasi-equilibrium relaxes 
to the equilibrium with the slower relaxation times r 2 . If T\ — r 2 , the intermediate relaxation step to the quasi- 
equilibrium disappears from (j4|), and it reduces to the Bhatnagar-Gross-Krook (BGK) model. 



B. Triangle entropy method 

For a practical implementation, explicit form of the functions f° q (M) and f*(M,N) are required. While for most 
of the cases, the equilibrium / z eq can be found explicitly either in a closed form or in a form of expansion, explicit 
construction of the quasi-equilibrium is case-dependent. Here we suggest a simple way to find quasi-equilibria in 
explicit form by perturbation around the equilibrium 0, Q. Let us assume that the equilibrium f eq (M) has been 
found explicitly, and that near the equilibrium 

f = n(M) + 6f*(M,N). (8) 

Expanding the entropy function H (fT]) at the equilibrium to quadratic terms, we obtain 

h {r + sf) = H(n) + j2 sfi + 1) + \ E jk + °( s f 3 )- ( Q ) 

The quadratic expansion of the entropy function f9| maintains convexity, and we find its minimizer 8f* subject to 
the linear constraints 

n 

J2 M k*Sf* = 0, 

2=1 

n 

J2 N lM = Nt-Np, (10) 

i=l 
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where 7V cq = N(f eq (M)) are the values of the non-conserved fields at equilibrium. Solution to this minimization 
problem has the form, 

(k M k N \ 

l + J2*sM si + J2xsN si , (11) 
s=l s=l / 

where Lagrange multipliers A s and \s are found upon substituting function (1111) into the constrains (|10p . and solving 
the (kM + few) x (kM + few) linear algebraic problem. Introducing a (fcjvf + £w)-dimensional vector of Lagrange 
multipliers, ( A, x) = ( Ai , ■ ■ ■ , A feM , Xi , . . . , Xfe„ ) , and matrices 

n 

(A MM ) kl = M ufrM u , k,l = l,...,k M 

i=l 
n 

(A M N)ki = j^Mkifj^Nu, k = l,...,k M , 1 = 1,..., fa 

(^)iW = ^Nkif^Nu, k,l = l,...,k N , (12) 



we find the solution in the matrix form, 



A A ( Amm Amn \ ( 



X J V A T MN A NN J \N- 



(13) 



where T denotes transposition, and N — N cq is the /cN-dimensional vector (Ni — iV^ q , . . . , Nk N — N^ q ). Note that 
the solution depends linearly on the deviation of the non-conserved fields N — N eq and non-linearly on the locally 
conserved M. For this quasi-linear quasi-equilibrium, the entropy production becomes 

Thus, with the use of the triangle entropy method, the kinetic model satisfies the entropy production inequality 
(both of the two quadratic forms in (|14|) are non-positive semi-definite) once the populations remain close to the local 
equilibrium. This is sufficient for most of the applications. 



III. QUASI-EQUILIBRIUM LATTICE BOLTZMANN METHOD 

We shall now derive a second-order time discretization for the generic kinetic equation Since the derivation 
only uses the consistency condition ([3]), it is equally applicable to exact quasi-equilibria and to those obtained by 
the triangle entropy method ([8]). Following [lfj, kinetic equations ((4]) are integrated in time from t to t + St along 
characteristics, and the time integral of the right hand side is evaluated by trapezoidal rule. Introducing a map 

fi^9i = fi-^Qi(f), (15) 

the result is written as 

g i (x + c i 6t,t + 6t)=g i (x,t)-LJ 1 \ gi (x,t)-f*(x,t)] - ^[f*(x,t) - f° q (x,t)}, (16) 

T2 

where 

26t 

Ul ~ 2n + 6t' 
f*(x,t) = f*(M(f),N(f)), 

f^(x,t) = fP(M(f)). (17) 

Note that the right hand side in (|16p still contains terms which depend on populations / rather than on the function 
g (|15p . In order to obtain a closed- form equation solely in terms of functions g, we first remark that, taking locally 
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conserved moments of the map (fT5|) we have (same as in the BGK case) M(f) = M(g), thus / cq (/) = f cq (g) in the 
second term in (fTB|) . For the quasi-conserved fields N, the situation is slightly different. Evaluating the moments N 
of the map (fT5|) . we obtain 

N(g) = N(f) + ^r(N(f) N^(f)). (18) 
Inverting this relation, and substituting it into (|16|). after simple transformations we obtain: 

gi (x + CiSt,t + St) = (l-wi)3, + a;i hjM /f(M) + wi ( I ^ IL ) fi(M,N% (19) 
where M = M(g), while AT' is evaluated according to the rule 

N '={ 1 -^) N ^ + 2^h- t Neq ^ ^ 

Equation (|19p is the basic second-order time stepping algorithm for the quasi-equilibrium lattice Boltzmann models, 
and is the main result of this paper. It is important to note the shift in the dependence of the quasi-equilibrium 
population /*, it does not depend just on N(g) but rather on a convex linear combination between N(g) and the 
equilibrium value N eq (g). Discretization in space depends on the problem at hand (in the simplest case, the lattice 
Boltzmann discretization is readily applicable). We shall now proceed with specific examples of QELBM. 

IV. EXAMPLES 

A. One-component fluid with a given Prandtl number 

Navier-Stokes equations for a one-component compressible fluid are characterized by the Prandtl number, Pr = 
(C p h)/k, where C p is the specific heat at constant pressure (C p — (D + 2)/2 is specific heat of ideal gas below), p 
is viscosity coefficient, and k is thermal conductivity. When a single relaxation time kinetic equation is used (for 
example, the BGK model), this results in a fixed Prandtl number. As we shall see it below, physical consistency 
of the QE models, the relaxation time hierarchy ([7]), implies that two models with two different quasi-equilibria are 
required in order to cover the entire range of Prandtl number. For the sake of concreteness, we consider a simple, 
weakly compressible flow model on the standard two-dimensional D2Q9 lattice [fl [l3| . The locally conserved fields 
M are density p, momentum density j , and pressure density p, 

n 

J2ih »i,«?}/< = {pj^p + p-'f}. (21) 

i=l 

To second order in the momentum, the equilibrium reads @, [l3j : 



9 2(1 -f) 



1 + ~T~ W [ ViaV ^ 2(T^|) — 5a \ 



(22) 



This model operates in a small temperature window T = p/p, around the reference temperature To — 1/3 (see Refs. 

[l3| for an estimate of the temperature window). We remark that, due to a low symmetry of the velocity set, the 
lattice BGK (LBGK) model with the equilibrium ((22)) gives PrLBGK = 4 [l3[ (not Pr = 1, as in the continuous kinetic 
theory or thermal lattice Boltzmann models with a higher symmetry 5, 8]). In the first case (Pr < PrLBGK), the 
specified slow variables N are the components of the heat flux q, 

n 

q a = ^ (via - u a ) (vj -u) fj. (23) 

4=1 

where u = j/p is the mean velocity. In the opposite case (Pr > PrLBGK), the slow variables A^ are the components 
of the stress tensor © 



2 

(v ia - U a ) (Vifi ~ Up) - —5 a() {v t - It) 2 



fi- (24) 
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Explicit expressions for the pertinent quasi-equilibria, f*(p,j,P,q) (case (|23)l . Pr < Ptlbgk) and f*(p,j,p,&) (case 
(fM| . Pr > PrLBGK) are easily found using the explicit formulas of the triangle entropy method, Eqs. (fTTj) and (fT5)) . and 
are not displayed here. Using the Chapman-Enskog method, we derive the hydrodynamic equations (Navier-Stokes- 
Fourier) for the density, momentum and temperature from the continuous QE models (UJ) with the corresponding quasi- 
cquilibria fi(p,j,p,q) and f*(p,j,p, 0) (note that the explicit form of /* is not required to perform this analytical 
computation, the consistency condition [3] suffices). The viscosity p and the thermal conductivity k coefficients thus 
obtained, imply the following Prandtl number: 

It = pT on , k = ~pT r 2 =► Pr = 4^ < Pr LBG K (25) 
2 t 2 

P = P T q t 2i k=\ P T t 1 => Pr = 4^ > Pr LBG K (El- (26) 

z n 

A comment to these formulas is in order: The hierarchy of relaxation times (n < t<i) implies a restriction on the 
range of admissible Prandtl number when a specified quasi-equilibrium is used. In the case where the heat flux is 
considered as the slow variable, the thermal conductivity (relaxation rate of the heat flux) is proportional to the slow 
relaxation time, k ~ T2, while p <~ n. In the opposite case when the stress tensor is chosen as a slow variable, the 
dependence is inverted, p ~ r 2 and k ~ ri . 

Implementation of the present QE models is a straightforward application of equation (|19|) . In Fig. [I] simulation 
of Couette flow between parallel plates at different temperatures is compared with the analytical solution. Diffusive 
boundary condition was used Agreement between simulation and analytical solution is excellent. 

2.5 I 1 1 1 1 1 1 1 1 1 1 
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FIG. 1: Steady state temperature profile in Couette flow between parallel plates moving with a relative velocity U, at the 
temperature difference ST. Lower curve: QELBM (gSJ with Pr = 0.710 (air); Middle curve: LBGK model with Pr = 4; Upper 
curve: QELBM (f24)l with Pr = 8. Symbol: simulation; Line: analytical solution. The Eckert number, Ec = (U 2 )/(C P 5T) is 
1.5. 



B. Binary mixture at a given Schmidt number 

Our second example is the isothermal binary mixture of ideal fluids, A and B, with particles masse toa,b- The locally 
conserved variables M are the densities of the components, Pa,b, and the momentum of the mixture j = Ja+Jbi where 
j A B are the momenta of the components. The efficiency of diffusion mixing versus viscous dissipation of momentum in 
various fluids is characterized by the Schmidt number, Sc = p/ (pDab), where p is the viscosity coefficient, Dab is the 
binary diffusion coefficient, and p = pa+Pb- Let us further introduce the molar fractions Xa,b — "a,b/("a+^b) where 
"a,b = Pa,b/toa,b is the number density of components, and the reference Schmidt number, Sc = mAB/(/^A^B) 
with TTiAB = (paPb) / (pa + Pb) the reduced mass density. We shall now consider a pair of quasi-equilibrium models 
with two different quasi-equilibria which cover the entire range of Schmidt number, Sc Sc. For concreteness, we 
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use the standard isothermal D2Q9 model. For the single-component case, the equilibrium f^ q (p,j) to second order 
in momentum is obtained by setting p/p = T = l/3in Eq. (|22[) . The equilibrium of the mixture is then described 
by the populations / A q Bl = f- q (pA,B,j)- 

The case Sc < Sc has been already considered recently [HI, [ill, and is mentioned here for the sake of completeness. In 
this case we consider the momenta of the components, j A B as the quasi-conserved variables N, and the corresponding 
quasi-equilibrium is immediately read off the equilibrium of the single-component fluid; to second order in momentum it 
* s /a Bi(PA,B> 3a b) = / 4 cq (/°A,B, 3 a b)- The continuous QE model recovers equations of the isothermal hydrodynamics 
that is, the advection-diffusion equations for the densities pa,b and the Navier-Stokes equation for the momentum j 
with jU = TiPo, £>ab = T 2Po > wnere Po — Tq(ua + it>b) i s the pressure at the reference temperature T . This 
model recovers Sc = (ti/t>j)Sc, and according to the hierarchy of the relaxation times ([7]), pertains to the fluids with 
Sc<Sc. 

In the opposite case, Sc > Sc, we choose the stress tensors of the components, ©a,b, as the quasi-conserved fields N 
(see Eq. I24j) . Corresponding quasi-equilibria, / A Bi (pA,B, 3, ©a,b), are constructed using the triangle entropy method 
in the same way as the case Pr > Ptlbgk- 

The lattice Boltzmann implementation is based on (TT9")) and an interpolation step, as explained in [n], [T§] ■ In Fig. 
[2] we present a simulation of diffusion of two fluids with a high mass ratio me = 500mA in a setup where initially a 
half-space is filled with the 90 — 10% mixture and the other half-space - with the 10 — 90% mixture. Agreement with 
the analytical solution is excellent. 




FIG. 2: Binary diffusion for the case Sc > Sc and ms/m/i = 500 at different time steps t. Symbol: simulation; Line: analytical 
solution. Triangles: t — 500; circles: t — 3000; squares: t = 6000; diamonds: t = 9000. 



V. CONCLUSION 

In conclusion, we suggested a systematic, physically transparent and realizable approach to constructing lattice 
Boltzmann models for hydrodynamic systems. All the models considered herein require only the choice of the quasi- 
equilibria appropriate to the physical context of the problem. Following the same pattern, it is straightforward to 
construct kinetic models for bulk viscosity and chemical reactions. 
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